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Abstract 

This paper discusses the source term treatment in the numerical solution of elliptic partial differential 
equations for an interior grid generation problem in generalized curvilinear coordinates. The geometry 
considered is that of a planar cross-section of a generic spiral-bevel gear tooth typical of a pinion in the 
OH-58 helicopter transmission. The source terms used are appropriate for an interior grid domain where all 
the boundaries are prescribed via a combination of Dirichlet and Neumann boundary conditions. 

New constraints based on the Green’s Theorem are derived which uniquely determine the coefficients in 
the source terms 1 . These constraints are designed for boundary clustered grids where gradients in physical 
quantities need to be resolved adequately. However, it is seen that the present formulation works satisfactorily 
for mild clustering also. Thus, a fully automated elliptic grid generation technique is made possible where 
there is no need for a parametric study of these parameters since the new relations fix these free parameters 
uniquely. 
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1 Introduction 

There has been a large amount of effort devoted to developing and enhancing the grid generation capability 
1 >2,3,4 through the solution of elliptic partial differential equations (pdes). The elliptic pdes used in the 
grid generation problems near boundaries are similar to the equations used in nuclear physics, diffusion- 
reaction problems, vortex problems, electric space charge problems, steady state heat transfer (conduction 
and convection) through long thin fins, etc. In the grid generation problems, these pdes contain appropriate 
source terms that control the distribution of grid points especially near the boundaries. In the literature, 
the elliptic pdes used for grid generation are erroneously referred to as Poisson equations which contain 
source terms that are functions of only the independent variables, whereas, in the pdes for grid generation, 
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these inhomogeneous terms also contain terms proportional to the dependent variables. Actually, in grid 
generation problems, close to a curvilinear boundary, the governing equations reduce to the long thin fin 
heat transfer equations with a finite heat transfer coefficient in the transverse direction (normal to the plane 
of paper) and a large heat transfer coefficient in the lateral direction. 

The focus in the studies referred to above has been on developing body conforming grids around bodies 
for external fluid flow simulations. The grids thus generated are smooth with at least first two derivatives 
continuous, appropriately stretched or clustered normal to any given coordinate direction and orthogonal 
over most of the grid domain. The inhomogeneous terms afford a grid control to satisfy clustering and 
orthogonality around specific surfaces (in three dimensions) and lines (in two dimensions). 

In external flows, these inhomogeneous terms, i.e., the source terms and the dependent variable propor- 
tional terms are designed to vanish away from the body so the problem reduces to solving a Laplacian away 
from the body. 

In the present study, the inhomogeneous terms used are appropriate for an interior grid generation 
problem where all the boundaries enveloping the grid will affect the solution through these terms. These 
terms are designed by interpreting their meaning physically through the principle of conservation of thermal 
energy close to the grid boundaries. 

The geometry treated here is that of a planar cross-section of a spiral-bevel pinion gear tooth typical of 
the OH-58 helicopter transmission pinion. This study is driven by the need to generate time-series vibration 
signatures from the OH-58 helicopter transmission by finite difference simulation of the appropriate structural 
dynamic equations. The choice of elliptic pdes for grid generation is entailed by the need to generate 
time series data as accurately as possible (see relative comparison with other representative grid generation 
methods in Ref. 5). 


2 Problem Definition 

The two-dimensional governing equations for an elliptic grid generation problem in an appropriately defined 
planar domain are 1,2 

+ iyy — P(£ y T}) 

T]xx +Vyy = QtiiV) 

where f and 77 are the generalized curvilinear coordinates, x and y are the Cartesian coordinates, and the 
P(£,rj) and Q(£,t)) are the inhomogeneous terms. 

The form of the inhomogeneous terms, P and Q, is, e.g., exponential 2 and is given by 

P(Z, *?) = “a*(77)«pn(C - ii)exp{-bi\Z - fc|) (la) 

Q{Z,V) = -Ci(£)sgn(ri - 7]i)exp(-di\T} - rfc|) (16) 

where i refers to the grid boundary in question. 

For the sake of argument, without loss of generality, if we take the case where f > & and 77 > 77* , then 
we have the inhomogeneous terms as 


p {& V) = -aifa)exp(-6i(£ - &)) (lc) 

Q(Z>V) ~ -Ci(Qexp{-di(rj - 77*)) (Id) 

At the boundaries, where £ = and 77 = 77^ Equations (lc) and (Id) respectively become 


= “<^(77) 
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<?(£,*?») = ~ c »(0 

When — &| or di\q — 77 , | is small, the inhomogeneous terms take the form given by 

P(t, V) = - &)) 


and 

Therefore, the governing equation for, e.g., £, in the vicinity of the boundary fi, becomes 

ixx + Zw = -a* (»/)(! “ “ €i)) 


ixz + Cyy - ai(r/)ft*4 = — - ai(q)b«& ( 2 ) 

If the term, a» 6 *£, were absent, the resulting equation would turn out to be a Poisson equation. The 
equation given above arises, e.g., in the steady state heat conduction problems in long thin fins, where £ 
is the temperature and where the heat transfer coefficient in the transverse thin direction is moderate but 
is large in the lateral direction, and there is a balance amongst the heat conducted through the fin, heat 
carried away from or to it through convection in proportion to this moderate heat transfer coefficient and 
the heat sources/sinks distributed over the domain. 

If we define a new variable 

0 = e - & 

then Equation (2) becomes 

6 XX + 0 VV “ a* W = -a* (3) 

The term, -a*, can be thought of as a heat source/sink term. 

Equation(3) tells us that when f > &, there is a balance between the heat convected from a control 
volume in the interior to the boundary heat conducted out of this control volume and the heat lost from 
the control volume due to the heat sink, <^( 77 ). Conversely, when f < &, there is a balance between the 
heat convected from the boundary £, to a control volume in the interior, heat conducted out of the control 
volume and the heat generated in the control volume due to the source, 01 ( 77 ). 

Prom Equation (3), it can be seen that for a given convective heat flux (given number of grid lines), as 
the product, aji>» decreases, the heat transfer coefficient decreases proportionally in magnitude which means 
that the temperature gradient at the boundary £* has increased so that £ approaches & rapidly. This means 
that there is a large gradient in £ from the grid boundary i to the interior, thereby resulting in a highly 
clustered grid near the boundary. 

Similarly, if we consider the case when & > £, then we have 

£xx + fw - a* (77)61$ = -OtO?) - G,i(ri)bi€i 
or, 

$ xx + 0yy - O,i(r))bi0 = ai{r]) 

where 0 = & - $ 

Away from this grid boundary, 6 *|$ — $*| or 6* (77 — T]i\ is large, and we are left with the Laplace equation, 
Af = 0 or Atj = 0 . Extremum principle is unconditionally maintained there, since the solution is harmonic 
in this case. 

Referring to Equation (3), the Green’s Theorem gives us 



where S is the surface area of a closed domain, C is the boundary enclosing this domain, n is the normal to 
the surface, da is the elemental area and ds is an elemental arc. 

The integrands on the left hand side, — a* and aib{9 represent the heat source/sink term and the convection 
term respectively, and the integrand on the right hand side represents the heat flux through the boundary 
0 . 

Equation(4) is used as a constraint to fix b t uniquely for a solution consistent with the specification of 
the boundary data. The extremum principle will be satisfied at the ith boundary, which is the requirement 
in the grid generation problems, since the energy conservation principle is satisfied. The term, is 

calculated iteratively through the solution process, which together with bi ensures the grid orthogonality and 
a given grid spacing at the ith grid boundary. 

In the design of these inhomogeneous terms, there is no restriction on the nature of the source term, 
di. It can change sign which indicates the presence of sources and sinks, subject to the constraint given by 
Equation(4). Otherwise, improper combinations of sources and sinks will violate the extremum principle. If 
over the domain, there is a net rate of heat generation due to the source/sink combination, then there has 
to be a positive heat flux convected away and vice-versa. This requirement will automatically be satisfied 
by Equation(4). 

If there is a point heat source present in the domain, the isothermals (temperature contour lines) will 
tend to cluster around it since the gradients in the vicinity of the source will be positive toward the source 
and high, depending upon the strength of the source, and conversely for a heat sink. Same argument applies 
to a line heat source and sink. By analogy, if the source term turns out to be positive over some parts of 
the domain, then the curvilinear coordinate lines will tend towards lines with higher coordinate values and 
vice-versa. 


3 Solution Procedure 

First, the boundary data are selected appropriate to the physics of the problem, so that the gradients in 
physical quantities can be resolved adequately. Since there is a symmetry plane and a rotational symmetry 
present in the present problem, the grid is reflected about this symmetry plane and then rotated around com- 
pletely about a moving axis of periodicity, thus substantially reducing the computational effort in generating 
the gear tooth grid. 

Then, by interchanging the independent and dependent variables, the governing equations to be solved 
in the computational space (£, rj) become 

ax u - 20x iv + 7 x„„ = - + <?(^,»j)x,) 

- 2 0y i „ + 7S/ W = ~J 2 ( P (£,»?)»* +Q(i,v)Vv) 

These equations are solved in the computational space using a line SOR relaxation algorithm where each 
coordinate line in one curvilinear coordinate direction is solved semi-implicitly using the Thomas algorithm 
for tri-diagonal systems. The inhomogeneous terms referred to above are designed and incorporated so that 
a desired grid behavior near the boundaries is achieved 6 . 

The inhomogeneous problem is solved using a technique similar to that of Ref. 3 by over-relaxing the 
inhomogeneous terms during the iteration process. The inhomogeneous terms used in Ref. 3 are well suited 
for external boundary value problems where they allow for clustering in only one curvilinear coordinate 
direction, normal to the body. But, in internal boundary value problems, inhomogeneous terms have to take 
account of the influence of the boundaries in both curvilinear coordinate directions. The inhomogeneous 
terms used here allow for clustering in both coordinate directions. 

The inhomogeneous terms, P(f,r?) and Q(£, 77 ), are evaluated at the boundaries in terms of the left 
hand side at each line relaxation sweep. Then outward from each boundary, the inhomogeneous terms are 
attenuated through an exponential function in each direction, as discussed above. In f direction, outward 
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from a given & boundary, this exponential term is of the form, — aiexp(— 6j|£ — &|), and, in 77 direction, it is 
of the form, -~aiexp(-bi\r) - r}i\). 

The boundary constraint given by Equation (4) is applied to a finite slender strip by evaluating the 
heat source/sink term and the convective flux term over the strip close to the boundary with the heat flux 
calculated around the strip. 


4 Results 

Figure 1 shows a finite-difference grid model of pinion and driven gears in mesh. The grids were generated 
automatically without any manual prescription of the free parameters. 



Figure 1: Finite difference grid model of pinion and driven gears in mesh 
Figure 2 below shows a close-up of the pinion gear and shaft. 



Figure 2: Close-up view of the pinion gear and shaft 
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The spiral-bevel nineteen tooth pinion gear typical of a OH58 helicopter is shown in Fig. 3 and the pinion 
shaft is shown in Fig. 4 below. 



Figure 3: Finite difference grid of the 19- tooth spiral-bevel pinion gear 



Figure 4: Finite difference grid of the pinion shaft 
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In Fig. 5 below, a cross-section of the pinion gear is shown. In Fig. 6, the corresponding cross-sectional 
grids for an individual pinion gear tooth and pinion shaft are shown. The tooth grid shows a desired clustering 
near and around the tooth. 

As is the case with generalized curvilinear coordinates, the cross-sectional grid for the shaft does not 
contain any geometric singularities. The polar coordinate singularity at the origin has been removed and 
manifests itself through the four triangular grid cells at the boundary as is seen in Fig. 6(b) below. 



Figure 5: Cross-section of the pinion gear 




(a) 


(b) 


Figure 6: Cross-sectional grids (a) pinion gear tooth (b) pinion shaft 
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5 Concluding Remarks 

The boundary constraints for elliptic grid generation problems developed in this study have been demon- 
strated to be applicable to a practical problem of gear teeth grid generation. Smooth clustered grids have 
been generated using these constraints without any recourse to redistribution of grid points which has been 
a common approach used in elliptic grid generation problems until now. With new constraints, elliptic grids 
can be generated in simulation time without any manual intervention thus making problems of structural 
dynamics and fluid dynamics over compliant boundaries straightforward. 
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